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The Dirac's chord method may be suitable in different areas of physics for the representation of 
certain six-dimensional integrals for a convex body using the probability density of the chord length 
distribution. For a homogeneous model with a nonconvex body inside a medium with identical 
properties an analogue of the Dirac's chord method may be obtained, if to use so-called generalized 
chord distribution. The function is defined as normalized second derivative of the autocorrelation 
function. For nonconvex bodies this second derivative may have negative values and could not be 
directly related with a probability density. An interpretation of such a function using alternating 
sums of probability densities is considered. Such quasi-probability distributions may be used for 
Monte Carlo calculations of some integrals for a single body of arbitrary shape and for systems with 
two or more objects and such applications are also discussed in this work. 
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I. INTRODUCTION 



Let us consider an integral 



> 

where Bi, B2 are three-dimensional bodies, ri — {xi,yi, zi) £ Bi, r2 — (2^2,2/2,22) £ B2 are pair of points, dVi — 
dxidyidzi and = dx2dy2dz2- 
'n] Similar integrals are used in different physical applications, e.g. in the calculations with a point-kernel method in 

, the radiation shielding and dosimetry P, Q • 

■ li Bi = B2 = B — is the single convex body, the Dirac 's chord method 3] may be applied for the calculation of the 
, particular case of the double integral Eq. ([1]) over pairs of points in the convex body B using the probability density 

■ of the chord length distribution, fi{l) 

^Biv)^ I I f l''^~''^l^ rfVidV2^^ r^ix)dxdr)dl. (2) 

K> , JbJb 'i^\ri-r2\'^ 4 Jq \JoJo ' 

^ . Here vectors ri, r2 £ B represent pair of points of the body B and Sb is the surface area of B. The infinite upper 
■ ■ ■ ' limit of integration is written for simplicity in the right-hand side of Eq. ([2]) and similar equations below due to 
obvious property fi{l) = for I > Z,nax, where /max is the maximal possible length of a chord. 

Such a formula may be used in analytical and numerical methods of the calculation of the integrals such as TfBi^)- 
A demonstrative advantage is the reduction of a six-dimensional integral to an easier expression such as Eq. ([2]). It 
is possible to obtain a direct analytical formula for the chord length distribution (CLD) for some bodies and it was 
initially used by Dirac et al, Ref. |j. 

Analytical expressions may be found only for few simple shapes and it is reasonable to consider application of 
Eq. ([2]) for numerical calculations of integrals, e.g. for Monte Carlo methods. Indeed, both the Monte Carlo method 
[3| and the Dirac's chord method from very beginning were used for the solution of analogue problems of the particle 
transport. The possibility to get rid of the singularity in the left-hand side of Eq. ^ is important for the 

application of Monte Carlo methods and it may be also actual for Eq. ^ with two neighboring or overlapping regions 
Bi and i?2- 
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However, even the generalization of Eq. ([2]) for a single nonconvex body is not obvious, because a straight line may 
intersect the body few times and an appropriate choice of a definition of CLD is not quite clear in such a case. There 
are three widely used nonequivalent constructions of CLD for a nonconvex body |6l-ll3l| . All intervals of the same line 
inside of a nonconvex body may be considered as separate chords to produce the multi-chord distribution (MCD). It 
is also possible to calculate the sum of lengths of all such intervals to define the one-chord distribution (OCD). 

The third definition introduces a generalized chord distribution as the second derivative of the autocorrelation 
function divided on some normalizer (e.g., Sb/^) [1, 113, lU- It is justified, because for a convex body such a 
formal expression is equal to the probability density for CLD. In a more general case such a definition is also useful, 
because just the generalized chord distribution should be used in Eq. ([2]) for a nonconvex body B instead of CLD 
[15'] and it is discussed below. However, for some nonconvex bodies the function may be negative for certain ranges 
of argument [O, [l^ . 

Methods of construction of such functions as alternating (in sign of terms) sums of probability densities are utilized 
in the presented paper. Such an approach provides a direct analogue of Eq. ^ for calculation of integrals for nonconvex 
bodies (15i] . An extension of this technique may be appropriate for treatment of a more difficult case Eq. ([T]) with two 
different bodies [iGj . 

Plan of the paper. In Sec. |TT] is revisited a ray method as a facilitated analogue of the Dirac's chord method. It 
produces an understanding physical model and introduces simplified versions of some tools applied further for chords. 
In Sec, mil some equations are collected which are useful further for discussion about applications of the chord method 
in Sec. lIVI The integral Eq. ([1} with two bodies and a multi-body case are discussed in Sec.|Vl Methods of applications 
of considered techniques for the statistical (Monte Carlo) sampling are discussed mainly in sections IIIBl IIVBI and 
IVSllVDl 



II. RAY METHOD 

A. Ray length distribution 

There is an analog of Eq. ([2]) with the probability density of the ray length distribution (RLD), l{1), i.e. instead of 
a full chord only a ray (segment) is considered. It is drawn from a point inside the body to the surface. The points 
have the uniform distribution and the directions of the rays are isotropic. It may be written in such a case 

Db{v)^Vb J^Lil)[Jjix)dx)dl, (3) 

where Vb is the volume of B. This expression could be considered as an intermediate step in the derivation of Eq. ^ 
in Ref. and might be simpler for explanation adduced below. 

Let us introduce a simple isotropic homogeneous model with particles emitted inside a convex body B and traveling 
along straight lines. If absorption of energy on the distance / from a source is defined by (p{l), the left-hand side of 
Eq. ([2]) or Eq. (jS]) with six-dimensional integral describes a fraction of energy absorbed inside the body. 

On the other hand, the same value may be calculated using the distribution of particle tracks (rays) inside the 
body. The part of energy, absorbed on a ray with a length I is 

Ivil) = / •f{x)dx (4) 
Jo 

and a fraction of rays with the length I is described by RLD l{1). It concludes an informal visual explanation of 
Eq. ([3]), because the total amount of emitted particles is proportional to the volume of B. 

The example with rays is also useful for the explanation of an appearance of alternating sums of distributions. Let 
us consider the nonconvex body and the ray with three intersections with the boundary depicted in Fig. [TJ 

For each such ray instead of Eq. (j4|) for the calculation of the energy absorbed inside the nonconvex body an 
expression 

nh rh 

/ ip{x)dx+ (p{x)dx = I^ih) - I^ih) + I^ik) (5) 
Jo J12 
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FIG. 1. Ray in nonconvex body 

should be used, where is the antiderivative of (p defined by Eq. It is possible to introduce few distributions 
tfc(Z) of distances from the source to fc-th intersection and to write instead of Eq. ^ 



fc=i -^0 Vo ^ 

Jo J Vo ^ 



(6) 



where fcmax is the maximal number of intersections of a ray with the boundary of B. The alternating sum in square 
brackets in Eq. ([6]) may be considered as a "quasi-probability distribution" 1(1) and so Eq. ([6]) may be rewritten to 
produce an analogue of Eq. ([3]) 



/.OO / V ^max 

^B{f) = Vb / m( / v{x)dx)dl, = ^(-i; 
Jo ^ 



fc+1 



^k{l). (7) 



More rigorous treatment may use so-called signed measures (charges) [l7| instead of term quasi-probability distribution 
used here. Some details may be found in Ref. 15. 

The visual interpretation of equations for rays above is rather informal. It was used understanding description with 
particles propagated along straight lines. Such a picture may create a wrong impression about impossibility to apply 
considered methods to more difficult models with scattering. It is not so, because the only essential condition is the 
possibility to use in integrals like Eq. ^ expressions depending merely on |ri — r2\. 

An example of appropriate model is a convex body with absence of scattering, but yet another case is an arbitrary 
body inside the medium with indistinguishable properties. The last case ensures possibility to apply Eqs. ([U [S]) and 
further generalizations to expressions with so-called build-up factors used in dosimetry and radiation shielding to take 
into account the scattering [ll, ^ . For the uniform and isotropic case such a build-up factor (for given energy) is again 
depending only on the distance from a point source. 

Let us introduce polar coordinates in the second integral in left-hand side of Eq. ([2]). It makes the consideration 
more rigour Q . Then for a convex body B it is possible to write 

^ dVi sindde / d(P / ^-^dR, (8) 

Jb Jo Jo Jo 47r 

where with the preceding notation of Eq. ^ R = \ri — together with 0,0 are polar coordinates of the vector 
R = r2 — ri and l{ri, 9, 4>) is the length of a ray from a point ri with a direction given by the polar angles 9 and cj). 
A designation dQ = sin9d9d(j) for the integration on the solid angle, 51 may be used for brevity 

S)b(v') = / '^Vi / — / ip{R)dR, (9) 

where Z(ri, r2) is the length of a ray from a point ri, i.e. denotes direction, represented earlier via 9, (j). 

Equation Q may be now derived, if to take into account normalizing multipliers Vb (volume of body B) and Att 
(area of surface of unit sphere). It is explained below in Sec. IIIII Some additional technical discussion and references 
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may be also found in Ref. [T^ Here is important to emphasize, that the ray in Eq. ^ is not necessary a particle 
trajectory, but a formal "axis" R of the integration on the variable R. 

Moreover, the formulas such as Eq. ([5|) or Eq. ([6]) with integrals on few disjoint intervals for a nonconvex body 
are also appropriate here and so Eq. ([7]) is valid. It justifies application of considered methods for arbitrary isotropic 
uniform media, i.e. for models with scattering. 

The important example is a body (convex or nonconvex) inside an environment with identical or similar properties. 
In such a case the term in the left-hand side of Eq. ^ depends only on distance |ri — even for points r2 near 
the boundary. For a convex body with straight tracks it is also true, but the environment does not matter, because 
trajectories of particles between two points inside the body may not fall outside the boundaries unlike the case with 
scattering. 



B. Method Monte Carlo with rays 

A useful application of Eqs. ([U [3]) is the Monte Carlo calculation of integrals. There is an additional advantage for 
the calculation of such integrals with many different (p{l) for each body. In such a case CLD or RLD for given body 
is calculated only once and used further with different functions fH). Functions, expressed via the definite integrals 
(single or double) of ip{l) in right-hand side of the equations may be calculated either numerically or analytically. 

The Monte Carlo sampling of a distribution is a standard procedure and may be visually represented as some 
histogram. The space between zero and the maximal possible length is divided on n bins, i.e. sections Ij < I < Ij -I- A^, 
J = 0, . . . , n — 1 and during simulation for each step an amount of "hits" in an appropriate bin is increased by one. 
For the equal size Al of all sections the index j of a bin is simply the integer part of l/Al and the tracing of such a 
data in the Monte Carlo simulations is fairly fast and useful procedure. 

For the apphcation of Eq. ([7|) it is possible instead of construction of k^^ax different distributions to create at 
once. If a ray intersects boundary in few points it is necessary to consider intervals from the origin to all points of 
intersection. For the length of each interval with odd index (first, third, etc.) it is necessary to add unit to number 
of hits in a bin, but for interval with even index it is necessary to subtract unit from a number in the relevant bin. 

Such a method describes the Monte Carlo algorithm for the generation of the function More difficult algorithms 
for quasi-probability distributions of chords are discussed below. However, it is reasonable at first to recollect some 
concepts for the explanation, why such algorithms are relevant with alternative definitions via derivatives of the 
autocorrelation function. 



III. HELPFUL ANALYTICAL EQUATIONS 



There are few functions related with presented models. It was already mentioned the chord length (distribution) 

density the ray length (distribution) density and the autocorrelation function, denoted further as 7(^). It 

is also convenient to consider the probability density of the distance distribution (DD) r]{l). There are important 
relations between these functions [6l-[l0l [H, [H, [H, [H, [l^ 

Kl) = -^7"(0, (10) 

MO = -^"^'(0 (11) 

^(0 = -f^7'(0, (12) 

2 



Vil) = ^7(0, (13) 

where Vb is the volume of body B and I = Jq°° I jjL(t)dl is the average chord length, that may be found for a convex 
body from a widely used Cauchy relationship [isl |20| - |22| 

r=4^. (14) 

OB 

The autocorrelation function ^{l) is defined here for body with density p{r) = 1 for r e _B as 

l{r) = J^pir,)p{r,+r)dV,, ^{l) = ^ ^^^^{r)dn, (15) 
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i.e. dVi = dxidyidzi, ri = {xi,yi,zi) and 7(Z) is the average of jjr) on a sphere with radius /. Definition of 7 here 
is lack of 1/Vb multiplier in comparison with some other works [la ] and it causes an insignificant difference in few 
equations. In fact, only the property Eq. (fT3|) is used further and the formal definition Eq. ([Tsl) is presented here for 
completeness. 

The probability density of the distance distribution r]{l) is easily defined for convex, nonconvex cases, and also for 
a system of two bodies. For the explanation of relations between derivatives of 7 in Eqs. ([TOl [T2|) it is convenient to 
start with the expression 

It may be explained using a statistical approach convenient here due to discussion on the Monte Carlo sampling. 
The left-hand side of Eq. ([TB| may be considered as an average of a function $(-R) = if{R) / [AttR^) of a variable 
R = \ri — defined on a space B x B. 

The multiplier is a measure (6D volume) of the six-dimensional space B x B and the division on this value is 
due to averaging. On the other hand, the standard correspondence [23j of a stochastic average with a mathematical 
expectation E makes possible to use the formula 

E$(i?) = j <^{l)diFR(l), (17) 

where Fail) is the cumulative distribution function of a random variable R and diFji{l) denotes probability density of 
R, but in the considered case it is just the density of the distances distribution (DD) defined earlier r]{l)dl = diF]i{l). 
Really, R is the distance between two points ri, r2 with independent uniform distributions and double integral in 
Eq. corresponds to averaging on B x B for such a pair. 
Equation (IT6l) may be rewritten due to relation Eq. ([T3| as 



Db(^)- / v{l)l{l)dl. (18) 
Jo 

After integration by parts it is possible to obtain from Eq. (IT5|) 

S)b(^) - /T-7'(0] ( fv{x)dx) dl. (19) 



For a convex body Eq. (|19p is in agreement with Eqs. ([51 [T^ . 

Equation ([3]) may also be proven using an analogue of the statistical approach discussed above. A detailed proof 



may be found elsewhere 15| and is only briefly sketched here. It is possible to consider Eq. (jH]) as an averaging on 
the five-dimensional space of rays, represented as product i? x S of the body B on the unit sphere S. It is necessary 
to use for normalization the volume Vb of B multiplied on 47r, the surface area of the unit sphere. In such a case 
Eq. ([3]) may be considered as an analogue of Eq. (jl7p for the mathematical expectation of some function depending 
on the length of a ray. 

It is possible to derive an equivalent of Eq. (jl2p for a nonconvex body with introduced in Eq. ([7]) if to use 
generalized functions and derivatives. The idea of a generalized function is convenient also for the further work with 
integrals such as S)((p). 

The generalized function (distribution) is defined [l7j as the continuous linear functional T((/)) on the space of test 
functions (f>. Usual integrable function may be associated with a functional defined for the test function (f>{x) as 



T^((/>) = / ^{x)<j){x)dx. (20) 
J —00 

On the other hand, "Db ■ ^ — ?■ "^si^) is also a linear functional on a test function ip and may be considered as some 
generalized function on (0, 00) defined by given body B. A topology on the space of test functions and continuity [l3l, 
i.e. DB{^k) 2)b(¥') for (fk f, are not discussed here. 

It is often used as a simplified notation tp instead of 1^ for such a regular generalized function Eq. (1^0]) p7| . In 
such a case Eq. (jl8l) may be rewritten simply as T> b = 1b- 

The generalized derivative [l7| is defined as functional 

1'icf,) = -1(0')- (21) 

Due to Eq. (jS]) it is possible for convex body B to write D'b — —Vb i-b and Eq. ([7]) for an arbitrary body ensures 
2)b = -Vb~lb- 
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IV. CHORD METHOD 
A. Chord length distribution 

For a convex body Eq. ^ may be rewritten with generalized functions and derivatives as 2)^ = {Sb/4:)^b- Here 
generalized functions may be more appropriate, because due to the expression with the second derivative CLD is not 
a regular function even if DD is not smooth in some points. 

Formally, for a convex body an expression with an additional integration along a chord appears due to a rearrange- 
ment of the integral Eq. ([8]) and consideration of all possible rays with origins along the same chord (see Fig.[|). 



FIG. 2. Consideration of all possible rays along a chord 

If to use a compact notation f dC for the formal integration on a four-dimensional space of straight lines [l^, [2l[ , 
it is possible to rewrite Eq. ^ after such a rearrangement as 

®b(<^) = / ^ / / 'f{x)dxdr, (22) 

where 1{C) is the length of a chord produced by the intersection of the straight line C with the convex body B. 
Equation is an analogue of a familiar expression used for the derivation of the Dirac's chord method (see Eq. (1.5) 
in Ref.i). 

In such a case there is a double integral along a chord due to the additional integration on sources of rays 

4'nO= / / V{x)dxdr. (23) 



Jo Jo 

For a nonconvex body and few chords, i.e. n intervals of intersections [x2kiX2k+i\, fc = 0, . . . n — 1 of the body by 
the same straight line, it is necessary to include only the integration on both "source" points r and "target" points 
r2 — r + X inside these intervals. Using rather technical calculation (see Ref. T^, Sec 3.3) it is possible to obtain 
instead of Eq. p3| the more difhcult expression 

I^\xo,...,X2n-l) = ^l!^\x2k+l -X2k) 
k=0 

n-1 k-1 

+ im[4'Ha^2fc+l - X2j) + l!^^\x2k - X2j + l)] 
k=l j=0 

- E E[4'^(^2fc+1 - X2, + l) + 4'Hx2k - X2,)]. (24) 
fc=l j=0 

It includes all n{2n — 1) possible ordered pairs Xk — Xj with indexes < j < k < 2n ~ 1 and may be rewritten as 

2n-l k-l 

l(J\xo,...,X2n-l)= E E(-l)'^'^'^f (^^-^^O- (25) 

k=l j=0 

E.g., for two intersections there are six terms (see Fig. [S]) 

ll'\xo, ...,X,)^ (Xl - Xo) + I^^\X, X2) 

+I!^HX3 - xo) + 4'HX2 - X^) - /(2)(X2 - xo) - 4'HX3 - Xl). 
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FIG. 3. Chords in nonconvex body and six possible segments 



Let us rewrite the integral Eq. ([22] 



(26) 



where for a convex body due to Eq. (|23)) l!^\C) = I^'{lc)- The same expression also may be used for a nonconvex 
body if to denote l!^'' (C) = I^'^ {xq, . . . , a^fn-i): where ccg , . . . , X2„_i designate all intersections of the straight line C 
with the boundary of B. 

(2) 

On the other hand, (£) may be expressed as a sum Eq. (^5]) with all possible (ordered) pairs of points. It is 
possible to rewrite Eq. ((26)) for the nonconvex case 



r(2) 



2ji-l fc-1 



(27) 



fc=l j=0 



The situation is similar with expressions for rays in a nonconvex body Eqs. (|6l[7]). Let us denote as ^J■jk{l) probability 
ensities of distril 
If to introduce 



densities of distributions of lengths Ijk = — produced by n{2n — 1) ordered pairs {xj" , x^) on a line C. 



where m is the normalization 



2n-l fc-1 

/^*o*(0 = E E(-l)'"'^V.fc(0, A* = m-i/itot(0, 

fc=l j=0 



iltot{l)dl, 



it is possible to write an analogue of Eq. ([2]) for a nonconvex body B 



^b{^) = SB ^ tlil)l!^\l)dl = Sb ip{x)dxdryi, 



(28) 



(29) 



(30) 



where is some constant. 

For a convex body fl{l) = s_b = 'S'b/4 and Eq. (|30l) may be explained using an idea with the averaging and 
the mathematical expectation Eq. (|T7)) already discussed for DD and RLD. Let us consider an average of the function 

(2) 

f{C) — 1^ (Ic) on the four-dimensional set C[B] of all straight lines intersecting the body B 



— f 4'Hic)dc^ ri!^\iui)di, 

WB Jc[B] Jo 



(31) 



where Wb is a measure (4D voluine) of C[B]. For a convex body it may be expressed as wb = t^Sb due to a 
Cauchy relationship [l^ [20l - [2^ [23 . and after comparison of Eq. (PT|) with Eq. (12^ we obtain necessary coefficient 
M;B/(47r) = 5*3/4 used in Eq. 

For a nonconvex body there are n{2n — 1) distributions fijkil) instead of one and Eq. pop is obtained via the 
alternating sum Eq. (j27p of these distributions and so formula sb = ttlb^wb / {4:Tr) is valid with wb is a measure for 
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a set of all straight lines intersecting B and is a constant used in definition of fl{l) Eq. ([28]) . The problem here 
is an absence of simple methods of a calculation wb and for nonconvex bodies and so it may be convenient to 
consider yet another approach for finding s b ■ 

It is possible to use an analogue of the relation in Eq. (ITT|) . The integration by parts of Eq. ([7]) for a nonconvex 
body produces 



SsM - Vb n-l'{l)][j j ^(x)dxdryi (32) 



f.1 

lo "-Jo Jo 

and after the comparison with Eq. ()30|1 it is possible to write 

-VBi'il) = SBflil) 

-Vb / ll'{l)dl = SB / ljl{l)dl 
Jo Jo 



Vb / lil)dl = SbI / fl{l)dl, 
Jo Jo 

where by definition Ib = J lfj-B(J)dl/ J jlB{l)dl and due to normalization of 1(1) and 

~SB = Vb/Ib- (33) 

For a convex body Ib is the average chord length. For a nonconvex body it is equal to the average chord length for 
the multi-chord distribution (MCD) mentioned earlier, because sums of lengths of all intervals in two last terms of 
Eq. ((24|) compensate each other 15]. It is clarified below in Sec. lIVBl about Monte Carlo simulations. 

In fact, the Cauchy relationship Eq. ([T4|) for the average chord length for MCD is proved for a broad class of 
nonconvex bodies [7] and so it is also possible to write due to Ec^. (1551) in such a case 

~sb - 5's/4. (34) 

For numerical methods using Eq. (j33]) with Ib sometimes may be preferable. It is instructive to consider an example 
with so-called voxel presentation of a body as a decomposition on small cubes or parallelepipeds. In such a case the 
surface is not smooth and a problem with correct approximation of surface area may not be resolved even for a formal 
limiting case with cubes of arbitrary small dimensions, e.g. for a sphere such a limit is Qirr^ instead of surface area 
47rr^. 



B. Method Monte Carlo with chords 



Let us consider some questions of the Monte Carlo generation of the quasi-probability distribution For each 

straight line with n > 1 intervals inside a body B it is necessary to consider 2n points of intersection with the boundary 
of B. Tangent points should be counted twice. Let us mark all points by real numbers Xfc, A; = 0, . . . , 2n — 1, there 
xq = and other Xk denote distances along a given straight line, i.e. Xk = l^^fc — '"o|, where denote positions of 2n 
points of intersections in the three-dimensional space. 

It is clear from further consideration that it is possible to use opposite order of points r2n-k-i and so 

"± orientation" of a line does not matter, i.e. two opposite directions could not be distinguished. Anyway, in applied 
calculations it is often more convenient to use directed lines. Standard algorithms of the generation of uniform isotropic 
(pseudo)random sequences of straight lines should be discussed elsewhere. 

Let us discuss the procedure of the construction of jl{l) for the given line. If the line intersects a body n times, 
it is necessary to consider set of 2n numbers Xk defined above and to calculate lengths Ijk = {xk — Xj) for all 
j,k:0<j<k<2n-l. 

For given Ijk a number in the relevant bin should be increased by unit for odd k — j and decreased by unit otherwise, 
i.e. if fc — J is even. For 2n indexes there are 1 -|- 2 -I- • • • -I- (2n — 1) = n(2n — 1) ordered pairs. Between them n 
pairs {x2k,X2k+i), k = 0, . . . , n — 1 represent usual chords lying completely inside the body and they have positive 
contributions. 

Remaining n{2n — 2) pairs are not forming intervals entirely belonging to the body and may be divided in two 
equal groups. There are n{n — 1) pairs with positive contribution, i.e. {x2j , X2k+i) or (a;2j+i, X2/c), 0<j<fc<n — 1. 
For other n{n — 1) pairs, i.e. {x2j,X2k) or {x2j+i, X2k+i), 0<j<fc<n — 1, numbers in appropriate bins should be 
decreased. 
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It is also clear from such a representation why sums of lengths of the pairs in two last groups compensate each 
other: 

{X2k+1 - X2j) + {X2k - X2j+l) = {X2k " X2j) + {x2k+l " X2j + l)- 

So sum of lengths alternating in signs is equal with summation of only n positive contributions due to chords of the 
straight line inside the body. 

There is also additional subtlety with normalization. For each straight line the total increase of values in all affected 
bins is 

^Ntot ^ n + n{n — 1) — n{n — 1) = n. 

So there are two different counters: the number of lines Ni and the sum of numbers in all bins Ntot ^ Ni^ which 
is equivalent with a total number of separate chords lying completely inside the body. The (quasi-probability) 
distribution should be divided in Nt^t for normalization. It is similar with the multi-chord distribution (MCD), 
because it is also normalized on the same number of separate chords Ntof. 

It was shown above that total sum of all lengths while taking into account signs is equal with sum of chords inside 
the body. But the normalization is the same as for MCD case and so the formal averaging of the variable I for the 
quasi-probability distribution fl{l) constructed here is the same as the average chord length for MCD that could be 
produced from the same set of straight lines. In a limit Ni ^ oo it ensures equality of / for both distributions already 
mentioned and used above in Eq. ([M)) . 

The total number of lines Ni corresponds to the normalization for OCD case, when for each straight line only one 
"aggregated" chord, equivalent to union of all chords inside a nonconvex body, is considered. Ni is also related with 
measure of set of all straight lines intersecting the considered body. Earlier in Eq. (j3ip this measure was denoted as 
wb- 

Yet another application of the both Ntot and Ni is the calculation of a constant m used earlier in Eq. ([^5)) . It may 
be expressed as a relation between fltot (that is not normalized on unit due to the contribution of lines intersecting 
the body more than one time) and ft. So m is the limit of ratio between total number of chords Ntot {cf MCD) and 
total number of straight lines Ni {cf OCD) 

m= lim Ntot/Ni. (35) 

Ni—^oo 



V. MULTI-BODY CASE 
A. Some equations with two different bodies 

This section is devoted to initial question about the calculation of the integral Eq. ([Ij with two different bodies. 
Here, it is also convenient to consider a simple model with particles moving along straight lines in isotropic uniform 
medium and to use the interpretation of Eq. ([T]) as a fraction of energy emitted in Bi and absorbed in B2 . 



•* *^ 




FIG. 4. Ray from Bi with interval inside of B2 

Let us consider a particle emitted in the first body with the straight trajectory intersecting the second one (Fig.U]). 
If the law of absorption is the same in both bodies and medium between them, it is possible to describe amount of 
energy absorbed in the second body as 

b 

ip{x)dx ^ I^{h) - I^{a), (36) 

where a and h are distances from a source to two intersections of second body by the ray and /(^ is defined above in 
Eq. Q. 
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Calculation of Eq. ([T]) for a convex B2 would be related with an analogue of Eq. (jS]) 

^B^M^I^dVil sin Ode I dcpj (37) 



where ^min, ^max, 0min, (/"max describe angular limits of integrations for given point ri and a, b are radial distances for 
given point and direction. It maybe simpler to use an analogue of Eq. (|9]) 

dY^ — ^{R)dR, (38) 



where S(ri, B2) is central projection from point ri of body B2 to surface of unit sphere. Yet, an alternative method 
of calculations is presented further and Eqs. (I37[ I38p are mentioned here rather for some clarification and comparison. 



B. Relation with methods for single body 



For two nonconvex bodies expressions above could be even more difficult, but it is possible to use a general principle 
to adapt already developed approach with single body [l6|. Let us choose both bodies as sources and consider four 
integrals ^g^{ip), s — 1,2, t — 1,2, i.e. — Ssi, — Ssai -SbI — 'SbI- Each integral takes into account only 
particles emitted in Bs and absorbed in Bt. 

The double integrals Eqs. ([1] [2]) comply with simple relations 

Sb.us. {^) = dll (^) + (^) + dll (^) + d'^l (^) (39) 

and 

2;?f J (ip) - Ds.uB, (^) - (^) - (^). (40) 

So many equations with two bodies may be reduced to already discussed case with the single body using union of 
these bodies B — Bi U B2. 

Here it is suggested that Bi does not intersect i?2. For overlapping bodies the decomposition on three parts: 
BiU B2, Bi\ B2 and B2 \ Bi should be taken into account. Instead of Eq. (|40p in such a case a modified expression 
may be used: 

2^?B? iv) = SsiUB. (<p) + ^B,nB, iv) - (v) - Sb. (<p). (41) 

Due to such equations, the computational methods discussed above make possible to find Eq. ([T]) after separate 
calculations of three or four terms in Eq. (|40l) or Eq. (HTj) . However, more direct approach discussed further is also 
useful and may be simply generalized for the case with many bodies. 

A simpler case of two disjoint bodies is suitable for almost straightforward modifications of Monte Carlo algorithms 
discussed above [T6| . Here Eq. p9p demonstrates that distributions obtained in simulation may be divided in four 
parts (for each pair source-target) without significant modification of algorithms for general nonconvex body discussed 
above and it may be even more convenient for explanation than Eq. (j40p . 

For two discontiguous bodies a union B = Bi U B2 formally is always nonconvex, because a line between a point 
in Bi and a point in B2 lies partially outside the union. So, here, consideration of quasi-probability distributions is 
especially justified. 

For non-overlapping bodies with adjoining boundaries B may be convex even with nonconvex Bi or B2. It corre- 
sponds to the case briefly mentioned below in a note about zones at the end of Sec. IV Dl It is enough to consider 
some body B and formally split it into two zones (parts) Bi and B2 ■ Even for convex body such parts may be either 
convex or nonconvex. These subtleties do not affect on methods discussed here. Insignificant change is necessary only 
for a case with overlapping bodies due to contribution of intersections with nonzero volume outlined in Eq. (I4ip . 



C. Application to calculations with rays 



Source points uniformly distributed in both bodies with equivalent density are used for sampling with rays. All 
intersections of rays with boundaries of both bodies are checked and appropriate bins are changed in four histograms 
Hst{l) marked by indexes source-target. 
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Joint consideration of all distributions lets to tackle a problem with normalization. For the function 1 term "quasi- 
probability distribution" could be justified due to normalization on unit integral and some relations with the proba- 
bility density for length of rays in a convex body, but if to write an analogue of Eq. ([3]) for the integral Eq. ([T|) 

SbI (^) = Wi2 J%2) (0 {fv{^)dx) dl, (42) 

where W12 is an unknown constant, it is simple to show that i(i2)(0 ma-y not be normalized for disjoint bodies, 
because due to Eq. pO)) 



W12 J l(i2){l)dl = {Vi + V2) J~iu{l)dl ~ Vi Jh{l)dl - V2 J~i2{l)dl = 0, 

where integrals of all functions t2(0 ^^'^ '■u(0 — '--BiUB2(0 ^-rc normalized on unit. 

However, if to include all four densities as components in the single process described by a (quasi-probability) 
distribution, introduced earlier 



{l) = J2hst)il), (43) 



it is possible to consider t(st)(0 as elements of some matrix 1(1) with the common normalization. The same approach 
may be used for more than two bodies. 



D. Application to calculations with chords 

The expression of Eq. ^ via chord distributions also may use similar principles [l^ . Here, it is also appropriate 
to use the decomposition Eq. (j39p . A straight line again determines 2n values xq,. . . ,X2n-ii but boundaries of both 
bodies must be marked appropriately. Each intersection should be refined as x^. with additional index s — \,2 for Bi , 
B2. 

Each pair (x|, x\) already has two additional indexes s and t representing four possible combinations for two bodies 

and it produces distributions combined with appropriate signs (— 1)*^+^^^ into jl[st){^)i s,t = 1, 2. It should be 

only mentioned that due to a symmetry for the chords "source" and "target" bodies could be hardly distinguished. 
Due to such property it is reasonable to use only three separate histograms H22 and H12 + H21 and to define 
M{st} (0 as a symmetric matrix. 

It may be directly generalized for a case with m bodies with s,t — 1, . . . , to. Advantages of application of discussed 
methods may be illustrated by consideration of a domain with many different bodies intersected by the variety of 
straight lines. It is possible to calculate all to^ integrals during the same Monte Carlo simulation. 

Here only to(to + l)/2 integrals are different due to symmetry, but anyway it may be a big number. For medical 
applications with 15 — 20 objects (organs) it is the calculation of hundreds values in a single Monte Carlo run. In 
fact, speed up may be even more critical due to possibility to split each body into few zones. The subdivision may 
be necessary for taking into account a variation of the intensity of emitters in the different parts of some objects. 

There is a subtlety for the calculation with few zones: it is necessary formally to split each boundary between two 
zones into two coinciding surfaces. In such a case all equations above are valid, but there are some intervals with zero 
length. Such intervals may be simply omitted, because integration along them produces zero values. 



E. Analytical expressions for two bodies 

There are useful analogues of expressions discussed in the Sec.|TTT]for the case with two bodies. A function r?(i2)(0 
may be defined as the probability density of distances between a pair of points uniformly distributed in first and second 
body, respectively. Analytical expressions written below may be used for the testing of the Monte Carlo simulation 
and some clarification. Technical details may be found in Ref. [l^ 

The correlation function 7(12) (0 defined for two bodies with unit densities Pfe(r) = 1 for r E Bk, fc = 1, 2 as 

7(i2)(r-) =^ Pi(n)p2(n+r)rfVi, 7(12) (0 = ^ ^7(12) (r-)df^. (44) 
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It is possible to derive the direct analogue of Eq. p3)) 

V(i2) = i^^J{i2){l) (45) 

and to write the generalization of Eq. ^TE\\ 

/ ^{i)jii2)m. (46) 

Jo 

Using integrations of Eq. (|46p by parts it is possible to express i(i2) and /i(i2) as first and second derivatives of 
correlation function 7(12) (0? respectively. It is similar with Eq. (1121) and Eq. (jlOl) . If to use the normalization on the 
union of bodies suggested above and the Cauchy relationship for the average chord length Eq. (|14l) , it may be written 

^(I2)(0 = ^^^|^7ri2)(0, (47) 

^a^)(0 = - ^^^|^^j (i2)(0. (48) 



Due to such normalization for the ray method an "unknown constant" in Eq. (|42p may be chosen as 

1^12 = Vb, + Vb, (49) 
and desired generalization of the Dirac's chord method for the initial equation Eq. ((T)) may be written finally as 

^B. ^ Sb, + Sb, n^^^ ^ tr^{x)dx dr) dl. (50) 



VI. CONCLUSION 



A novel approach with the application of quasi-probability distributions (signed measures) to calculations of integrals 
such as Eqs. ([T][2]) is advocated in presented work. It may be useful in many areas of physics. This paper is written with 
a purpose to present a fairly brief, but closed description of considered methods. Additional technical details, proofs 
of some equations together with appropriate links with theory of geometrical probabilities may be found elsewhere 

[il[i3. 

It is shown, how models with ray and chord length distributions suitable for a single convex body should be altered 
for a nonconvex case and multi-body systems. An essential new property of such extensions is the necessity to use 
instead of probability densities some functions, which sometimes do not satisfy the non-negativity condition. 

Maybe such a counterintuitive "negative probability" produced certain difficulties and a delay in development 
and applications of these methods despite of high effectiveness of numerical algorithms based on ray and chord 
distributions. On the other hand, qua si-p robability distributions are rather common in quantum physics after so- 
called Wigner function representation [26| and Feynman wrote an essay about the concept of negative probability 
with reasonable examples both in quantum and classical physics p7j . 

In fact, the functions and do not necessarily directly related with probability distributions and so should 
not cause some conceptual challenges. Appearance of negative values may be simply illustrated using Eq. ([5]) and 
Fig. [H Here the ray (0,^3) includes a ray (0,Zi) already taken into account and the interval (Zi,/2) outside of the 
body, that should not be counted at all. 

For the work with an interval {I2, 13) expressions such as Eq. ([5|) were used, but it may be described in the standard 
probability theory. If probability measures are known for sets Ri — A, R2 = AU B , R3 — AU B U C , it is possible to 
write for C: P(C) = P(i?3 \ R2) = P(i?3) - P(i?2) and for AuC: P{A U C) = P(i?i) - P(i?2) + P(i?3). 

Overlapping sets, i.e. rays with the same origin, are used in construction of (see Fig.[T]). Positive and negative 
terms such as P(i?3) and — P(i?2), used for the calculation of the same P(C), affect two ranges of argument 1(1). So 
for I = i?3 there is some positive gain, but for I = R2 there is corresponding decrease and it may produce negative 
values of t{l) for some intervals of /. An extra hit is added to some bin. A removal from another bin — is an effort 
to compensate that, but it may produce a negative result. 

The construction of is simpler, than the generalization of the chord length distribution /i(Z), but a reason of 
appearance of negative values in both cases is similar. An amount of terms in expressions for a ray grows linearly 
with respect to a number of intersections and for chord it is quadratic dependence. The structure of sets is also more 
complicated for chords, but here alternating signs in formulas such as Eq. (j24p again correspond to an expression with 
unions and differences of some overlapped sets. 
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